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ESTIMATION  OF  TIME-OF-ARRIVAL 
OF  UNDERWATER  ACOUSTIC  SIGNALS 
BY  SPLINE  FUNCTIONS 


INTRODUCTION 

In  the  study  of  underwater  acoustic  signals,  one  of  the  most  important  problems  is 
to  determine  the  timc-of-arrival  (or  signal  onset)  from  the  given  data.  For  instance,  the 
transducer  and  panel  transfer  function  identification  techniques  all  depend  on  an  accu¬ 
rate  estimate  of  this  value.  In  this  report,  the  underwater  acoustic  signal  will  be  repre¬ 
sented  by  a  spline  curve  whose  initial  knot  lies  in  the  interior  of  the  time  interval  so  that 
the  spline  function  is  identically  zero  to  the  left  of  this  knot  and  “takes  off'  to  the  right 
at  this  knot.  Hence,  the  initial  knot  clearly  defines  the  time-of-arrival  of  the  acoustic  sig¬ 
nal.  The  objectives  of  this  research  are  to  study  a  mathematical  model  in  the  form  of  an 
extremal  problem  whose  solution  yields  an  accurate  estimate  of  the  time-of-arrival  and 
to  develop  an  algorithm  for  solving  this  extremal  problem. 

The  first  section  will  be  devoted  to  the  study  of  spline  functions  with  special  em¬ 
phasis  on  the  computational  procedure  of  each  polynomial  piece  exactly  and  determining 
the  coefficient  matrix  in  a  modified  (or  “penalized”)  least-squares  problem.  Both  the 
treatment  and  results  here  are  different  from  those  in  the  literature.  When  the  discrete 
least-squares  problem  is  solved,  the  coefficient  matrix  is  usually  singular.  We  will  choose 
the  solution  whose  least-squares  measurement  is  minimized.  This  minimum  solution  is 
achieved  by  taking  the  so-called  Moore-Penrose  pseudoinverse  of  the  coefficient  matrix,  a 
topic  that  will  be  discussed  in  the  section  entitled  The  Minimum- Normed  Leas,  Squares 
Solution.  This  presentation  is  intended  to  be  complete  so  that  both  known  an<  new  re¬ 
sults  are  included.  The  mathematical  model  that  determines  the  initial  knot,  or  equiva¬ 
lently  the  time-of-arrival  of  the  acoustic  signal,  will  be  posed  and  studied  in  the  section 
entitled  Estimation  of  Tiine-of-Arrival.  An  algorithm  to  determine  the  solution  of  the 
corresponding  extremal  problem  is  described.  A  computer  program  with  numerical  ex¬ 
amples  will  also  be  included. 

FITTING  OF  UNDERWATER  ACOUSTIC  SIGNALS  BY  SPLINE  CURVES 

Underwater  acoustic  signals,  noisy  or  not,  can  best  be  fitted  by  using  spline  curves. 
In  addition  to  the  usual  benefits  from  spline  representation  such  as  efficiency  in  compu¬ 
tation,  flexibility  in  choosing  the  order  of  smoothness  vs  minimizing  the  possibility  of 
oscillatory  behaviour,  the  variation  diminishing  property,  etc.,  the  main  reason  in  us¬ 
ing  spline  curves  to  fit  underwater  acoustic  signals  is  that  the  initial  knot  of  the  knot 
sequence  of  the  approximating  spline  function  clearly  defines  the  time-of-arrival  (or  sig¬ 
nal  onset).  This  section  is  devoted  to  the  study  of  spline  functions.  Although  there  is  a 
vast  literature  on  this  subject,  we  will  introduce  a  new  approach  which  cannot  be  found 
in  any  published  paper  or  book,  except  partially  in  Refs.  3  and  5.  in  order  to  facilitate 
our  self-contained  discussion  and  to  improve  the  computational  procedure  to  suit  our 
purposes.  In  particular,  since  the  initial  knot  is  going  to  be  an  interior  knot,  the  results 
in  this  section  are  somewhat  different  from  those  in  the  spline  literature.  This  section 
is  divided  into  five  subsections.  Since  a  spline  function  is  a  piecewise  polynomial  func¬ 
tion,  the  first  subsection  contains  a  short  discussion  of  representation  of  a  polynomial  by 
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‘‘Bernstein  coefficients.”  One  advantage  is,  of  course,  that  a  Bernstein  polynomial  de¬ 
fined  by  using  the  more  general  formulation  here,  is  coordinate  independent,  so  that  it 
provides  an  ideal  representation  of  any  polynomial  piece  of  a  spline  function  between  two 
consecutive  knots.  Spline  functions  are  defined  in  the  second  subsection  with  empha¬ 
sis  on  the  ones  with  uniform  knot  sequences,  and  a  computational  method  introduced 
in  Ref.  5  using  the  Bernstein  coefficients  directly  is  included  in  the  subsection  entitled 
Computation  of  £?-Splines.  The  last  two  subsections  are  important  for  the  remaining 
material  of  this  report.  The  new  idea  of  representing  an  underwater  acoustic  signal  by 
a  spline  curve  with  the  first  initial  knot  defining  the  time-of- arrival  is  introduced  in  the 
subsection  entitled  Spline  Representation  of  Underwater  Acoustic  Signals.  A  modified 
(or  penalized)  least-squares  procedure  is  used  where  the  modification  is  achieved  by  in¬ 
troducing  a  parameter  that  governs  the  noise  level.  The  final  subsection  includes  a  com¬ 
putational  method  of  the  coefficient  matrices.  This  information  is  especially  important 
for  our  computational  task. 

Representation  of  Polynomials 

Let  k  be  a  non-negative  integer  and  tta-  denote  the  vector  space  of  all  real  polynomi¬ 
als  of  degree  at  most  k.  That  is,  each  p  £  tr*  is  of  the  form 

k 

P(x)  =  J2aJxl  (1) 

j= o 

where  <70,  •  ■  • ,  Ok  are  real  numbers.  Of  course,  the  collection  { 1,  x, ... ,  xk }  of  monomials 
provides  a  basis  for  zr* .  However,  for  both  theoretical  and  computational  purposes,  it  is 
usually  necessary  to  focus  our  attention  on  a  certain  interval  [a,  ft],  and  in  doing  so,  the 
basis  {l,.r, . . .  ,.Tk'}  is  no  longer  useful,  since  it  does  not  reflect  upon  this  interval.  For 
this  reason,  we  will  use  the  so-called  barycentric  coordinate  of  the  interval  [«,  b]  defined 
as  follows.  Let  u  and  v  be  two  linear  polynomials  defined  by 


u  :■=  u(x)  - 


b  —  x 


v  :=  e(.r)  = 


b  —  a 


Note  that  the  pair  (»,  v)  of  “variables”  identify  the  interval  [a,  b\  in  the  sense  that  u  rep¬ 
resents  the  initial  end-point,  a,  relative  to  the  interval  [a.  b],  and  v  the  final  end-point,  b. 
relative  to  this  interval,  namely: 


u(a)  =  1 
u(b)  =  0 


and 


v(n)  =  0 
v(b)  =  1. 


This  pair  of  variables  can  be  used  in  place  of  the  variable  x.  since 


x  =  ua  -f  vb.  ( 3 ) 

Of  course,  we  must  remind  ourselves  that  u  and  v  are  related  by  the  identity,  u  +  r  =  1 
for  all  x.  Now,  it  follows  that  the  collection  of  k  -f  1  polynomials 


where  0  <  i.j  <  k  and  ?  +  j  =  k ,  is  also  a  basis  of  that  is.  every  polynomial  ]>  d  ~ 
can  be  uniquely  represented  as 
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P(u,v)  =  p(u(x),u(z))  =  aijiPij(u^v)  (5) 

i+j—k 

for  all  x,  where  the  summation  is  taken  over  all  i  and  j  with  the  restrictions: 

0  <  i,j  <  k  and  i  +  j  =  k.  In  other  words,  the  coefficients  },  i  +  j  =  k.  called 

the  Bernstein  coefficients  of  p(u,v)  uniquely  determine  the  polynomial.  In  Fig.  1,  we 
display  the  Bernstein  coefficients  of  three  cubic  polynomials  on  [a,  6],  which  are 

/b  —  x\ 3  /x  -  a\ 3 
\b  —  a)  6  —  a) 

n/b  —  x\  (x  —  a\2 

3(— Kn)  • 

and 


(b~x) 

(2  ( x  — 

1  -  3( 

f b~x) 

(x~a\ 

\b- a) 

'  U-aJ 

\b  —  aJ 

\b-aJ 

respectively.  It  is  important  to  remark  that  not  only  the  collection  { a ktJ  } ,  i  +  j  =  k, 
uniquely  determines  p(u,  v),  but  it  also  gives  a  geometric  description  of  the  curve  traced 
by  the  polynomial  p(u,v).  Indeed,  since  >  0  for  all  x  G  [a,  6]  (or  equivalently,  0  < 
u,r  <  1)  and  =  1,  p(u,v)  is  a  convex  combination  of  its  Bernstein  coefficients  a 
and  hence,  lies  in  the  “convex  hull”  of  the  set 


<6> 

This  two-dimensional  set  is  called  the  Bernstein  net  of  p(u,r^.  By  joining  this  net  by 
straight  line  segments,  we  note  that  the  graph  of  p(a,v)  lies  in  the  “convex  hull”  as 
shown  in  Fig.  2,  where  both  the  Bernstein  nets  and  the  polynomial  curves  of  the  exam¬ 
ples  given  in  Fig.  1  are  shown. 


1  0  0  1 


a  b 


0  0  1  0 


a  b 


0  1-10 


a  b 


Fig.  1  -  .^-Coefficients  of  cubic  polynomials 
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Spline  Functions 

Let  k  be  a  non-negative  integer  and 

t:  —  <  *_i  <-•-<#,<  — 

be  a  bi-infinite  sequence  with  tt  — >  oo  and  — >  -x  as  ;  — >  oo.  Suppose  that  c  —  to 

and  cl  =  when  n  is  a  positive  integer.  Then  a  function  /  E  d]  (i.e..  /  has 

continuous  derivatives  up  to  order  k  —  1  on  the  interval  [c.  d})  is  called  a  spline,  function,  of 
order  (k  +  l)  on  [c, </]  with  knot  sequence  t,  if  the  restriction  of  /  on  each  interval 
is  a  polynomial  of  degree  at  most  k,i  =  •  ■  • ,  —  1. 0. 1.  ■  ■  We  will  denote  the  space  of 
spline  functions  of  order  (k  +  1)  on  fc, d]  with  knot  sequence  t  by  Sk,t(c,d). 

It  is  well  known  that  a  basis  of  the  space  Sk,t(c.d)  is  given  by  the  collection  of  so- 
called  5-splines  5*.  t  where  i  =  —  k,... , n  —  1,  a  totality  of  n  -f  k  functions.  Each  5- 
spliue  Bk  t,,  is  in  Sk, t{r,d)  and  vanishes  identically  outside  the  interval  (tl,t,+k+ 1  )■  The 
interval  [tt,t,+k+i J  is  called  the  support  of  the  5-spline  Bk.tj ■  it  is  also  well  known  that 
the  support  [/,,  fj+t  +  i]  's  minimal  in  the  sense  that  any  function  /  in  Sk ,t(c,  d)  whose 
support  is  a  proper  subset  of  [t,,  t{+k+i  ]  must  be  the  zero  function  [2],  In  Fig.  3.  we  give 
the  graphs  of  5-splines  of  order  1,  2,  3,  and  4. 


Fig.  3  -  5-Splines  of  order  1-4 

In  this  report,  we  are  only  concerned  with  spline  spaces  having  uniform  knot  se¬ 
quences  t;  that  is,  we  only  consider  t,+i  —  t,  =  t,  —  f,_i  for  all  i.  The  construction  of 
5-splines  with  uniform  knot  sequences  is  particularly  easy.  We  start  with  the  special 
case  where  t,  =  i. 

Let  N0(.r)  be  the  characteristic  function  of  the  interval  (0.1);  that  is. 


fl  if  0<c<l 
1  0  otherwise 


(7  a ) 


We  define,  inductively. 


=  /  AVi(-r  -  t)dt. 

Jo 

k  =  1.2 .  The  graphs  of  N0 .  N\ .  .V2  and  A't  are  shown  in  Fig.  4. 


( 7/)  i 
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Nq  N]  N'j  N3 


Fig.  4  -  5-Splines  at  integer  knots 

It  is  easy  to  show  that  TV*.  £  Ck~x{  —  00,  00),  and  since  each  integration  increases  the 
degree  of  the  polynomial  pieces  by  one,  we  see  that  Nk  is  in  Sk,z ,  where  Z  is  the  set  of 
all  integers.  The  B-splines  in  Sk,z  are  given  by 

Bk,z,i(x)  =  Nk{x  -  i)-  (8) 

Now  suppose  that  —t,  =  h>  0  for  all  i.  Then  the  B-splines  in  Sk.t^ic.d).  where 
the  knot  sequence  is 

with  /0  =  c  and  tn  =  d,  are  clearly  given  by 

Bk,tH.i(x)  =  Nk(l(x-c)-i).  (9) 

These  are  the  B-splines  that  will  be  used  in  this  report. 

Computation  of  B-Splines 

Let  t,+  j  —  ti  =  h  >  0  for  all  i  and 

t/,:  •••  <  t,  <  •••  <  ti  <  ■■■ 

with  to  =  c  and  tn  =  d.  In  view  of  the  formula  in  Eq.  (9),  to  obtain  the  B-spline  basis 
{Bk,th.i}'i  —  —k,. . .  ,n  —  1,  of  the  spline  space  SfctA(c,d),  it  is  sufficient  to  determine 
Nk(x).  Although  Nk ( x )  can  be  computed  by  using  the  definition  in  Eqs.  (7a)  and  (7b), 
it  is  more  efficient  to  compute  the  Bernstein  coefficients  of  each  polynomial  piece  by  fol¬ 
lowing  the  method  described  in  Ref.  3.  Since  N\(x)  is  piecewise  linear  with  values 

*■<‘>={0  ‘f  r,1.  iu,) 

where  i  €  Z  (Fig.  4),  its  Bernstein  coefficients  are  {0.  1}  and  {1,0}  on  the  intervals  [0.1] 
and  [1,2],  respectively,  as  shown  in  Fig.  5. 


0  1  0 


Fig.  5  -  B-Spline  A’i(.r) 
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To  determine  the  Bernstein  coefficients  of  AN(.f).  we  need  an  intermediate  step 
where  the  Bernstein  coefficients  of  ( l/2)[A"i  ( -r )  —  Xj(x  —  1)]  are  recorded.  Now  along  the 
interval  [0.1],  first  input  the  zero  initial  condition.  Then  add  this  zero  value  to  the  first 
value  of  (l/2)[iV1(x)  —  AN(J  —  1)],  namely  zero,  to  get  the  second  zero  value.  Next  add 
this  zero  value  to  the  second  Bernstein  coefficient  of  ( 1/2)[AN (x )  —  AN (.r  —  1 )].  namely  1/2. 
to  get  the  1/2  as  the  third  Bernstein  coefficient  of  AN(-r)  on  [0.1].  For  the  interval  [1.2j. 


we  perform  the  same  operation:  add  the  initial  value  of  1/2  to  the  first  value,  namely 
1/2,  for  ( l/2)[i'V1(.r )  —  N\(x  —  11],  to  get  the  second  Bernstein  coefficient  of  AN(.r)  on  [1. 
which  is  (1/2)  +  (1/2)  =  1.  ana  next,  add  this  value  1  to  the  second  value,  namely  —1/ 
for  ( 1  /2)LV]  (.r )  —  AN  ( x  —  1)],  to  get  the  third  Bernstein  coefficient  of  AN ( .r )  on  [1.2]  which 
is  1  —  (1/2)  =  1/2.  For  the  interval  [2.3].  we  repeat  the  same  procedure.  ( See  Fig.  0). 


Fig.  6  -  Construction  of  AN ( .r ) 

To  compute  AN ( .r ) ,  we  again  need  the  intermediate  step  of  (l/3)[AN(.r)  —  AN ( x  -  1 )] - 
We  then  follow  the  same  procedure  as  above  in  the  computation  of  AN (•'»')■  but  this  time 
the  support  of  AN ( x )  is  [0.4],  one  unit  longer.  ( See  Fig.  7). 


Fig.  7  -  Construction  of  AN(-t') 

In  Fig.  8,  we  show  the  computational  procedure  for  the  quartic  and  quintic  splines 
X\(x)  and  AN (x),  and  for  simplicity,  we  have  omitted  the  arrows  and  write  down  the 
Bernstein  nets  for  AN(.r)  -  AN(t  —  1)  and  AN(.r)  —  AN(Jr  —  1 ).  instead  of  (l/4)[AN(.r)  — 
AN(.r  -1)]  and  (l/5)[AN(-c)  —  AN(.r  —  1)],  respectively. 
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Fig.  8  -  Construction  of  JV4(:r)  and  iV5(x) 

Spline  Representation  of  Underwater  Acoustic  Signals 

The  spline  model  will  be  used  in  this  report  to  represent  underwater  acoustic  sig¬ 
nals.  Let  the  time  interval  in  signal  measurement  be  [0,d],  If  t0  >  0  is  the  time-of-arrival 
of  an  acoustic  signal,  then  the  spline  curve  that  represents  this  signal  must  be  identically 
zero  for  t  <  t0  and  “takes  off’  at  t  —  to.  Hence,  the  spline  curve  is  given  by  the  spline 
function 


n  —  1 


2—0 


(11) 


where  the  knot  sequence  is 


1ft'  to  \  <C  ■  '  '  tn  —  i  *C.  •  •  *  <C  tn-ffc 
with  0  <  t0  =  c,  tn  =  d,  and 


t>  =  tj-i  +  h  ,  h  =  - - ,  (12) 

n 

j  =  1,2,  ...,n  +  k.  (See  Fig.  9.) 

Here,  due  to  the  nature  of  the  spline  series  S*(t)  in  Eq.  (11),  this  spline  model  as¬ 
sumes  zero  values  of  5^(t0),. . . ,  1  ( t0 )  at  the  time-ot  arrival  t0..  Hence,  if  the  acoustic 

signal  should  have  nonzero  slope  at  the  time-of-arrival.  only  the  linear  spline  model  ( k  — 
1)  can  be  used.  In  the  forth  coming  report,  we  will  allow  stacked  knots  of  spline  func¬ 
tions  of  degree  k  at  t0  in  order  to  make  use  of  higher  degree  spline  models  even  though 
the  slope  at  the  time-of-arrival  may  not  be  zero.  In  doing  so,  the  procedure  becomes 
adaptive  in  nature  but  involves  more  complicated  computations. 
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Fig.  9  -  Spline  model  of  signal 

Let  \is  assume  that  the  measurement  is  taken  at  { r,},  where  0  <  rx  <■■■  <  t\  <  d, 
and  the  signal  measurement  is 


f(Ti)  =  fi  ■  i  = 

In  practice,  we  have  N  n.  The  coefficients  {e^}  in  the  spline  series  Sf.(t)  in  Eq.  (11) 
must  be  so  chosen  that  the  quantity 


11/ -  St||2  +  A  (13> 


is  minimized,  where  ||  ||  is  a  norm  to  be  specified,  and  A  a  non-negative  parameter  that 

is  selected  to  adjust  the  "smoothness"  of  the  acoustic  spline  curve  S^(t).  For  noise-free 
signal  {  f,}.  A  should  be  chosen  to  he  zero,  hut  a  positive  A  value  is  requied  if  the  signal 
is  contaminated  with  noise.  In  other  words,  the  value  of  A  must  he  adjusted  adaptively 
according  to  the  noise  level.  A  detailed  study  will  he  given  in  the  next  report  [4], 

Two  practical  norms  |(  ||  are  the  I.1  and  f2(w)  norms,  defined  respectively  by: 


Ml/.*  :  = 


MO  \*(lt 


1/2 


(14) 


and 


I !  fi  1 1  / 5 1  w  i  • 


( 15) 


where  W  =  { (/’, }  is  a  given  sequence  of  positive  numbers,  called  wenrlit^.  If  the  signal  is 
noise-free,  or  if  the  variance  of  the  noise  process  is  time  invariant,  we  may  take  ir,  --  1. 
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On  the  other  hand,  for  noisy  data  with  time- varying  variance,  {u!,j  should  be  chosen 
according  to  the  standard  deviation  of  the  noise.  When  the  L 2  norm  is  used,  the  signal 
function  f(t.)  must  be  measured  for  all  values  of  t  £  [0.  d\.  This  is  usually  not  feasible. 

and  one  way  to  get  around  it  is  to  define  f(t )  as  a  piecewise  linear  function  such  that 

/( r, )  =  /,,  i  =  1 . N ,  and  that  the  restriction  of  f(t)  to  each  interval  [r,.r,+1]  is  a 

linear  function.  1  <  i  <  N. 

Let 

,</  "-i 

F(c)~Fh(  c)=  / 

J0  ,_n 


For  the  time  being  let  us  assume  that  h  is  fixed  (which  is  equivalent  to  saying  that  the 
time-of-arrival  is  assumed  to  be  known).  Here, 


To  determine 


Sk(t)  =  J2cJBkM,J(t)  (1 

j=0 

which  minimizes  the  quantity  in  Eq.  (13),  or  equivalently  F(c)  in  Eq.  (16),  we  simply 
determine  the  following  “ normal  equations"  by  differentiating  F(c)  with  respect  to  c: 


su 


)ii  +  X£'-Jo 


f(t)BkM,i(t)dt, 


i  =  0, . . . ,  n  —  1.  Hence,  by  setting 


A-k ,h  =  f  -Bt-,th,j(f)-0*,tA,i(O^ 

.  -  0<t,j<n  —  1 


Jo  f(t)Bk.th.oindt 


L  Jo  f(i)Bh.th.n-\(t)<lf  J 

the  system  of  normal  equations  Eq.  (19)  can  be  written  in  the  matrix  form: 

(  Ak,h  +  J'In  )C  =  f h  ■ 

where  /„  is  the  order  identity  matrix  and 


9 


Chui 


L'O 


Cn  __  j 


(23) 


One  advantage  of  this  approach  is  that  since  Ak  k  is  the  Gramian  matrix  of  the 
B-splines  and  A  >  0,  the  coefficient  matrix  (Ak  k  +  A/n)  is  nonsingular,  so  that  c.  and 
hence  the  solution  Sk{t),  can  he  uniquely  determined. 

A  more  standard  approach  is  to  use  the  f2( w)  norm,  so  that  the  discrete  data 

{( r, ,/,)},  z  =  1, . .  .  ,  N,  can  be  input  directly  without  introducing  /.  The  main  advan¬ 
tage  of  this  approach  is  that  noisy  data  can  he  treated  much  more  easily  especially  in 
the  case  when  N  n.  Let 


N 


n  —  1 


n  —  1 


(7(c)  :=  Gh(  C)  =  y;  /,  -  CjDkithj(Tt)  wt  +  A  y  c) .  (24) 

1=1  j= 0  j  =  () 

Again,  for  the  time  being  let  us  assume  that  h  is  fixed.  Then,  to  determine 

n  —  1 

*,.>«)  (25) 

J  =  0 

which  minimizes  the  quantity  in  Eq.  (13),  or  equ  ently  G'(c)  in  Eq.  (24).  we  follow 
the  same  procedure  as  above  to  derive  the  following  normal  equations: 


n  —  I  /  S 


N 


y:  [y.  BkMtJ{T()BkM  t{Tt)wt  j  c* + Ac  •  =  y,  f'D *  t„.i(77)uy, 


f=i 


i  =  0, . . . ,  n  —  1;  or  in  matrix  form: 

T  A/„]c  —  f hi 

where  Ak ,h  is  the  n  x  n  matrix  whose  (i,j  )th  entry  is  given  by 


N 


y  BkXh,(T()Bk'th 


f=i 


and  the  vectors  c*  and  f/,  are 


n-i  J 


and 


t  /,  = 


.frBk.th.o{T(  )irf 


T.L  f'D*  .th  n  -  1  (Tf)ir, 


(2G) 


(27) 


[2S: 


(29) 
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respectively.  The  disadvantage  of  this  approach  is  that  the  coefficient  matrix  (A^.h  +  A/n ) 
is  frequently  singular  so  that  the  normal  equation  Eq.  (27)  might  not  have  a  unique  so¬ 
lution.  The  standard  approach  to  determine  a  solution  c*  of  Eq.  (27)  is  to  consider  an¬ 
other  least-squares  problem: 


min  IK-4*,/,  +  AJn)c  -  f*||#2. 

C 

However,  even  this  problem  usually  does  not  have  a  unique  solution.  Hence,  we  will  con¬ 
sider  another  extremal  problem,  by  choosing  c*  to  be  the  solution  of  the  above  least- 
squares  problem  with  the  minimum  ||  ||/>2  value.  It  turns  out  that  this  “minimal  so¬ 

lution’’  is  now  unique  and  can  be  determined  by  finding  the  so-called  Moore-Penrose 
pseudoinverse  of  (Akth  +  A/„).  This  topic  will  be  discussed  in  the  section  entitled  The 
Minimum- Normed  Least  Squares  Solution. 


Computation  of  the  Coefficient  Matrices 


To  solve  the  normal  equations,  Eqs.  (22)  and  (27),  it  is  necessary  to  determine  the 
corresponding  coefficient  matrices.  Hence,  the  quantities 


(30) 


and 


N 

*ij  =  X  BkMti(Te)Bk<th,j(Te)w(  (31) 

fei 

must  be  computed.  In  this  report,  to  produce  a  computational  efficient  algorithm,  we 
have  chosen  the  knot  sequence  t  to  satisfy  Eq.  (12),  so  that  the  material  presented  in 
the  subsection  entitled  Computation  of  R-Splines  can  be  applied.  In  particular,  by  Eq. 
(26)  we  have 


Bk^)  =  Nl‘{j(t  ~c)  -  *) 

=  Nk(n-i  -  ^{d-  «)), 

so  that 


(32) 


b>j 


Nk(t)Nk(t  +  i  —  j)dt\ 


ij  =  0, . . .  ,n  -  1. 


(33) 


where  we  have  used  the  property  that  ;V*(f)  =  0  for  /  <  0.  Now  the  restriction  of  .Y k(t) 
on  each  interval  [r,  r  -f  1],  where  r  =  0, 1 . . . ,  is  a  polynomial  of  degree  k  whose  Bernstein 
coefficients  can  be  computed  by  using  the  procedure  outlined  in  the  same  subjection. 
Hence,  to  determine  ,  we  are  led  to  compute 


where 
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and 


Pk(u,v)=  Y  c(m  («**'•) 

(+m  =  k 


(34) 


Qk(u,v)=  Y  dl,vtg(u’v)  (35) 

p  +  q=k 

are  two  (Bernstein)  polynomials  of  degree  <  k  on  [r,  r  +  1].  The  computation  is  straight¬ 
forward,  since 


_  C  A-!  T(f  +  p  +  l)r(m  -f  Cjf  -f  1) 

Clrn\  p\q\  F(C  +  p  +  w+q  +  2) 
kl  k\  (t  +  p)\(m  +  q)'. 
f!m!  p!<7?  (f  -f  m  +  p  +  q  +  1 )! 

_  (^!)2  (f  +  />)!(m  +  g)\ 

(2A:  +  1)!  (\m\p\q\ 

Therefore,  we  have 


y'  v-  (f  +  p)!(u;  +  q)'  k  k 

(2k  -hi)!  ^  ^  [\m\p\q\ 

e+m  —  k  p-f  ?=/t  *  * 


(37) 


where  c*m  and  djf/  are  the  Bernstein  coefficients  of  Pk  and  Qk  respectively. 

For  example,  if  a  linear  spline  curve  is  used,  then  by  using  the  Bernstein  coefficients 
of  N i ( .r )  in  Fig.  5  and  formulas  in  Eqs.  (33)  and  (37).  we  have 


-ft 

3 

for 

i  =  j  -  n  -  1 

9 

-ft 

3 

for 

*  =  i  =  0... 

. .  n  -  2 

I/, 

6 

for 

1*  -  =  1 

0 

otherwise. 

(3S) 


If  a  cubic  spline  curve  is  used  to  fit  the  underwater  acoustic  signal,  then  the  Bernstein 
coefficients  of  A \(.r)  in  Fig.  7  and  formulas  in  Eos.  (33)  and  1 37 j  can  be  used  to  com¬ 
pute  b3tj.  0  <?._/  <  n.  as  follows: 

(a)  For  i  =  j  and  i  =  Q, ...  ,n  —  4.  we  have 


;  »  ___  ■s-'  ((  +  p)\(6  -  (  -  p)'.  , 

"  ~  ~h 

where  (c„ . c3  )  =  (o,0,Q.  1/g)  and  (d„ . r/-,)  =  (l/G.  2/G.  4/G.  4/g)  . 


or 
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2  (3!)2  r  6! 


6!  2!4!  4!2! 


l3  2  (3! )  f  6!  (-  z!4!  A  i  4!2!  ^  b!  ^ 

6,1  =  62  7717  l  (3TP  +  L(3!)2  +  (2!)24  +  (2^2 16  +  (TTp1 


5!  2!4!  (3!)2 

4.  0 - 0  4  9 - 4  4  2- — —4 

2!3!“  3!2!  (3!)2 


(3!)2  4!2!  5!  4  2416, 

+  2(2iy  8  +  22!3!S  +  22!3!  1G1  i  h  =  7!  ^ 

(b)  For  i  =  j  —  n  ~  3,  we  have 


i-3 

3,n  —  3 


(3!)2  ^  v''  +  P)-(6  —  i  —  p)'r  ,  ,  ,  u 

—  E  E  Tijs-ows-p)!^'  +  2<W'!'' 


(c)  For  2  =  j  =71  —  2,  we  have 


; ,3  _  V'  +p)!(6  -  i  -  p)\ ,  ,  ,  1 , 

n_2,n_2  7!  ' pl 

1208, 


(d)  For  i  =  j  =  n  —  1,  we  have 


fe3  W  (7  +  P)!(6  ~  £  ~  P)!  r 

fe„-2,n-2  7!  1.  Z.  £1(3  -  £)!p!(3  -  p)! 1  pJ 


(e)  For  i  =  0, . . . ,  n  —  4  and  j  =  i  +  1,  we  have 


>3  w|2  il  +  p)!(6  ~  £  ~p)7n  j  ,  j  ill 

0  7!  EE£!(3_()!p!(3_p).l2'  f  "1  ' 


/here  ( d0,dud2,d3 )  =  ^4/6, 4/6, 2/6, 1/6^;  or 


,3  =l(3!)!f2f(3l)!  +  ^2+  14+A4 

•2  62  7!  t  l(3!)2  3!2!  3!2!  (3!)2 

6!  ,  2!4!  0  4!2!  „  6! 

4  - 44 - 84 - 8H - 4 

L(3!)2  (2!)2  (2!)2  (3!)2 

+  3^!12+i|18+(lt?17 

(3!)2 _ 4!2>  5'  n  . 


4!Z  0  'l 

20+M71S+^12]}'' 
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(f)  For  2  =  n  —  3  and  j  =  ti  -  2,  we  have 


63  -  V  V 

&n-3.„-2  -  -=r 


(f +  />)!(C  -  f  -P)! 

w-nw-p)' 

1062 


+  <l((lp]h 


7! 


-  /? . 


(g)  For  ?’  =  ri  —  2  and  j  =  n  —  1,  we  have 


i,3  _  (3!)2  y-'  ((  +  />)!(6  —  f  -  p)\ 

n~2'n~x  71 

129 


[cfdp}h 


7! 


/2. 


(h)  For  2  -0 . it  4  and  j  =  2  +  2.  we  have 


6?..  = 


(3!) 


((  +  P)'-(6  —  l  —  p)! 


->■ J  V-  \r-  i~  P):(  o  - 

7!  7!{3 -/)!/>! 


it;w-ow3-p)! 


l<(dv]h 


1 

—  — 0 


(3!)\  412!.  o: 


7! 
120 


7! 


5!  G1 

4  +  - 4  +  - o  4. 

1  OlOl  ’  olr.1  “  I 


(3!)2  2!3!  2!3!"  (3!)2 


fi)  For  2  —  22  —  3  and  j  =  22  —  1.  we  have 


fc 


3 

n  —  3 ,  n  —  1 


(3! )2  y"  y'  ^  +  —  />)! .  .  . 

7!  Sip;  2!(3-Oy(3-,.)! 


(j )  For  2  —  0. ....  22  —  4  and  j  =  2  +  3,  we  have 


h 


3 

<} 


'3! )2  ((  +  /-»)!( 6  —  (  —  /i)! 

7'  s,r,,!(3-ow3-p)! 

(3!)2  ( 3! )( 3! )  1  1 

7!  (3!)(3!)  6  2  '  7! 


where  f  I'd  .....  r-;{ )  =  (l/6.  0.  0.0  V 

fk)  For  j  2  —  j  |  >  4  .h'j  =  0;  and  b')t  =  by  the  definition  in  F<j.  i  30 ) . 

Summarizing  the  above  results,  we  may  write  down  the  matriees  .4^  h  fOI  /• 
3  ns  follows: 
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'4l  ft  ~  6 


4  1 

1  2. 


Al  i  —  37 


where  C„_3  is  an  (7*  —  3)  x  (n  —  3)  banded  Toeplitz  symmetric  matrix.  D  an  (77  —  3)  x  3 
matrix  with  transpose  Dr ,  and  £3  a  3  x  3  matrix  given  by: 


-2416 

991 

120 

1 

991 

2416 

991 

• 

0 

120 

991 

• 

•  * 

« 

1 

• 

• 

• 

• 

1 

C'„-3  — 

• 

« 

% 

• 

120 

• 

« 

*  * 

991 

0 

• 

• 

• 

1 

120  991 

2416 

120  1 


L  991 

120 

lJ 

2396 

1062 

60 

1062 

1208 

129 

60 

129 

20 

To  compute  6*-,  we  note  that  again  bk-t  ~  bk}  for  all  i,j  =  0, ...  ,7?  —  1  and 
<v 

i>ij  =  Y  Nk  (n  ”  *  ”  fr(d  -  Tt))Nk  (w  -  j  -  j{d  ~  r()j  (4*1 

(=\ 

where  again  A\( x )  can  be  computed  by  using  the  procedure  outlined  in  the  subsection 
entitled  Computation  of  B-Splines.  Note  that  .Y*(.r)  vanishes  outside  the  interval 
[0,  A-  4-  1],  so  that  most  of  the  terms  in  the  summation  in  Eq.  (44)  are  zero.  In  partic¬ 
ular.  we  again  give  the  Bernstein  coefficients  of  the  linear  B-spline  X\(.r)  and  the  cubic 
B-spline  iV3(.r)  in  the  following: 


Xi(x] 


0  0  0 
0 


1  1  1  1  1 


0  0  0 

T» 


V,.  r 


Chili 


THE  MINI  MUM- NORM  ED  LEAST  SQUARES  SOLUTION 

When  discrete  ch.ta  information  is  given,  a  standard  norm  to  use  for  the  quantity 
described  in  Eq.  (13)  to  be  minimized  is  the  (2(  w)  norm.  For  w  =  { ir, }  with  ir,  =  1.  we 
simply  write  ll  =  /2(w).  Hence,  for  any  sequence  c  =  {e*},  we  have 


As  shown  in  the  subsection  entitled  Spline  Representation  of  Underwater  Acoustic  Sig¬ 
nals,  the  corresponding  extremal  problem  then  reduces  to  a  linear  system  described  by 
Eq.  (27).  where  the  coefficient  matrix  (A*./,  -f  A /„)  is  frequently  singular.  Hence,  the  lin¬ 
ear  system,  Eq.  (27),  does  not  have  a  unique  solution  and  is  even  “numerically  inconsis¬ 
tent”.  To  overcome  this,  the  usual  method  is  to  minimize  the  t2  norm  of  the  difference 
(Ak  h  +  A/„  )c  —  f For  convenience,  we  will  simplify  the  notation  by  setting 

A  =  A*,/,  +  b  =  f/, .  (4G) 

Hence,  we  will  consider  the  extremal  problem: 

min  ((Ac  —  b jj /-2 .  (47) 

C 

Of  course,  if  the  original  system  Ac  =  b  is  consistent,  then  the  minimum  value  in  Eq. 
(47)  is  zero,  and  a  solution  to  the  extremal  problem  in  Eq.  (47)  also  solves  the  linear 
system  Ac  =  b  as  required.  In  any  cast',  whether  the  linear  system  Ac  =  b  is  consistent 
or  not,  there  is  no  guarantee  of  a  unique  solution  to  the  problem  in  Eq.  (47).  For  vari¬ 
ous  reasons  such  as  stability  (when  n  is  very  large),  the  desirable  solution  to  Eq.  (47)  is 
one  whose  f2  norm  is  also  minimized.  That  is.  we  will  consider  the  problem: 

min{||c||p:  |(Ac  -  b||/2  =  min  f| Ac  -  b||,s }.  (48) 

C 

We  will  then  choose  the  solution  of  Eq.  (48)  to  be  the  solution  c*  in  Eq.  (27). 

In  this  section,  we  will  see  that 

c*  =  A+  b  or  c*  =  ( AkJl  +  A /„  )+  f* .  (49) 

whe  re  A  +  is  the  so-called  Moore- Penrose  pseudoinver.se  of  A. 

Singular  Value  Decomposition  and  Moore-Penrose  Pseudoinverse 

The  singular  value  decomposition  of  an  arbitrary  matrix  is  studied  in  this  section. 
We  will  take  a  somewhat,  unusual  route  in  introducing  this  familiar  concept  so  that  the 
definition  of  the  Moore-Penrose  pseudoinverse  becomes  very  natural.  We  will  also  prow 
that  this  definition  is  independent  of  the  (non-unique)  singular  value  decomposition 
Let  .4  be  an  m  x  n  matrix  which  may  not  even  be  a  square  matrix  (although  in 
our  application,  it  is  always  square  and  symmetiic).  Suppose  that  rank  .1  k,  whole 
k  <  mini w,n).  Then  A1  A  and  AA1  are  n  x  n  and  in  x  in  (respectively)  non  ne^att'. 
definite  symmetric  matrices  of  rank  k  and  ha\e  the  same  eigenvalues 

n\  >  •  •  •  >  n[  ;.»  ()  ■■•-(). 
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which  we  arrange  in  non-increasing  order.  Let 

be  the  eigenvectors  of  ,4rd  corresponding  to  the  eigenvalues  {a\, . . .  ,a2k.  0, ....()}.  so 
chosen  that  they  form  an  orthonormal  set.  In  other  words,  the  n  x  n  matrix 


V'  =  [vj  . .  .  v„] 


is  unitary,  and 


a 


2 

1 


ArAV  = 


o 

ak  y 

0 


Next,  let  a,  be  the  positive  square  root  of  af  and  set 


Then  we  have 


or  equivalently, 


In  addition,  we  have 


u,  =  — A\t  .  i  =  1, . . . ,  k. 


At  Uj  =  —  (AtAv,)  =  —  a] 

<7,  (7j 


v,  =  —ATiii  .  i  =  l,...,fc. 


AAt  u,  =  aiA\l  =  a]  Uj, 


(50) 


(51 ) 


(52) 


(53) 


(54) 


so  that  u,  is  an  eigenvector  of  AAT  corresponding  to  the  eigenvalue  o\,  where 
i  =  1 _ _  k.  Let  Ujt+i , . . . ,  um  be  orthonormal  eigenvectors  corresponding  to  the  zero 

eigenvalues  of  AAT\  that  is, 


and 


T 

u  u, 


'•J 


.4/1tu,  =0,  t  =  A-  +  l,...,m, 
i.j  =  k  A  1, _ w.  Here.  6tJ  is  the  Kronecker  delta: 


(  1  for  i=j 
\  0  for  i  f  j. 


Now.  since  .4 .4 ^  is  symmetric,  eigenvectors  belonging  to  distinct  eigenvalues  are  orthogo¬ 
nal.  so  that 


u/  u7  =  0  for  i  =  1 . k  and  j  =  k  +  1 . in. 
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Next,  for  i,j  =  1 . k,  we  also  have 


1 


ata} 


'!  (a!~Avj)  =  -i-v.Vjv, 

\  /  (7,(7 , 


=  h 


since  (vi,.  .  .  .  vjt }  is  an  orthonormal  set.  Summarizing  the  above  properties  of  { n i 
u ,n } ,  we  see  that  the  m  x  m  matrix 


U  =  [uj  ur 


is  unitary  and 


aatu 


(7t 


0 


cr, 


o 


u. 


We  also  have  the  following. 

Lemma  1.  The  matrix  A  has  the  following  “ singular  value  decomposition" : 


(55) 


(5G) 


.4  =  UZVT 


(57) 


where 


£  = 


0\ 


o 


o 


Ok 


o 


m  x  n 


0  :  0. 

Proof.  S  nice  V  is  unitary,  Eq.  (57)  is  equivalent  to  .41'  =  IW.  or  equivalently 

.4v,  =  a,U] . 4v/t  =  akuk 

since  .4v*  +  1  =  •  •  •  =  .4vn  =  0  and 


(5S) 


\  O 


9) 


CE  =  [(7 1  u i  .  .  .rrkUk  0.  .  .()]. 

This  completes  the  proof  of  the  lemma  since  Eq.  (59)  is  the  same  a.'  Eq.  (51). 

In  view  of  the  result  in  Eq.  (57)  and  the  fact  that  l’1  —  V  ~  1  and  Y}  —  l'~'.  we 
are  now  ready  to  give  a  definition  of  the  "inverse"  of  the  matrix  .4  >  which  is  not  neces¬ 
sarily  square  and  not  necessarily  of  full  rank). 
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Definition.  The  Moore-Pcnrose  pseudoinverse  of  an  m  x  n  matrix  .4  is  given  bv 

A+  =  VE+Ur  (60) 

where 

o 

ak' 

o 

Since  the  unitary  matices  U  and  V”  are  not  unique,  we  must  prove  that  the  defini¬ 
tion  in  Eq.  (60)  is  independent  of  the  choices  of  U  and  V .  We  have  the  following. 

Theorem  1.  Let  .4  =  UTV1  =  U E V  r  be  two  singular  value  decompositions  of  A:  that 
is, 

-4.4 1  u,  =  of  u,  ,  -4-4ru,  =  of  Uj  ,  i  =  1, . . . ,  m 

AtAv,  —  ofv ,  ,  AtAv,  =  ofv i  ,  i  =  1 . n 

and 

u/uj  =  ujxij  =  6ij  , 

V,7  VJ  =  v^Vj  =  <5t;  ,  ij  =  1, . . . ,  n, 

where  17  =  [uj  . . .  um],  U  =  [uj  . . .  um],  V  =  [v,  . . .  v„],  and  V  =  [v,  . .  .  v„].  Then 

-4+  =  VE+UT  =  VZ+UT 
where  £+  is  defined  in  Eq.  (61). 


Proof.  Write  crj cr*  =  Aj, . . . ,  Aj, . . . ,  Af, . . . ,  A*,  where  Ai  >  >  Af  and 

n  u 

i ]  +  •  •  •  +  i(  =  k.  For  each  j  =  let  Pj  be  an  i}  x  i}  unitary  matrix  such  that 

— t-  —  i  “t- 1  •  •  •  h  j i  -f — f  ij )  i  -t — f  •  j  - 1  + 1  •  •  •  u ,  j  -j — *- 1;  \  Pj 

In  addition,  let  P(+ .  be  another  (n  —  k)  x  (n  —  k)  unitary  matrix  such  that 


[fifc+l  .  .  .  fi/TJ  -  [ujt-f  1  .  .  •  Unj  ]  Pf  ]  . 


Then  we  have 


V  =  u 
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Since  by  Eq.  (53),  we  have: 


1  iT- 

VP  =  TA  UP  •  P  =  A+-- 

•  +  ij- 1  +  1 

'  i  T  •  ■  ■  +  ij. 

. (.  it  then  follows  that 

,  1  ,7 

V,H - h'J-i-hl'--',iH 1-  *j  J  ^ 

- y 

[uu  +  ■  +tJ  _  i 

+  I  •  • 

■  «»o.  4  ••  +,,! 

+  1  •  - 

■  Hi,  -i - hi  j\Pj 

=  (VM  +  ' 

1  +  1  '  ‘  ' 

vn  + 

-+.,1  r, 

again  by  Eq.  (53). 

Let  Q  be  an  (n  —  k)  x  (n  —  k)  unitary  matrix  such  that 

[v*+i  ...v„]  =  [v*+1  . ..  v„]Q. 

Hence,  we  have 


This  yields 


(VZ+CT)r  =  UZ+VT 


or  1  as  required.  Here,  we  have  used  the  fact  that  PjPj  —  Itj.  This 

completes  the  proof  of  the  theorem. 

Characterization  of  the  Moore-Penrose  Pseudoinverse 

The  Moore-Penrose  pseudoinverse  ,4+  of  an  in  x  u  matrix  A  can  be  computed  by 

finding  the  eigenvalues  of  A1  A.  the  corresponding  orthogonal  eigenvectors  {rj . e„  } 

and  the  null  space  of  AA1  since  the  other  eigenvec  >rs  { #/  j . ft  a  }  of  A  A1  can  be  com 

puted  by  using  Eq.  (51  ).  In  the  following,  we  give  a  charartc  li/atioii  of  ,i+  in  terms  of 
some  matrix  identities. 
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Theorem  2.  .4+  is  the  Moore-Pcnrosc  psncdainverse  of  A  if  and  onlv  if  A+  satisfies  the 
following  properties: 

(i)  AA+A  =  A, 

(ii)  <AA+)r  =  44+ , 

(iii)  4+44+  =  4+,  and 

(iv)  (4+4)r  =  4+4. 

Proof.  Ir  is  clear  from  definition  that  .4+  satisfies  the  four  properties  listed.  Now  sup¬ 
pose  that  B  satisfies  (i)  through  (iv),  and  we  have  to  prove  that  B  =  ,4+  To  do  so 
recall  that 


where 


,4  =  U 


Si  0 
0  0 


Si  = 


<7\ 


o 


V‘ 


0 


We  define 


i?9  B\2 

Hi!  B22 


=  vtbu. 


Then  by  using  the  property  (i)  for  B,  we  have 

ABA  =  A , 

so  that 


(UTAV)(VTBU)(UTAV)  =  UtAV 


S,  O' 

Bg 

■s,  O' 

's,  o' 

.0  OJ 

B2 1  B22 

,0  0. 

.0  0. 

f  S,Z?9S,  =S, 

1  ^5,2=0  . 

Therefore,  Z?9  =  and  BV2  =  0. 

By  using  the  property  (ii)  for  B.  namely:  {AB)r  =  4B  or  (UrABU)r  =  f'7  4RT 
so  that 

[(r/'.4V')(r7'JBf-)]7'-(r/4r,(r/z?f-, 
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s.  oi  [sr1  o 

0  0  B2 1  Dn 


O  O  I  B'2  i  B 


we  nave 


'h 

SiB-Ji' 

h  o' 

0 

O 

B2lSx  0. 

so  that  Z?2i  -  0.  Similarly,  by  using  property  fiii)  for  D.  we  also  have  D> 2  =  0.  Henri', 
we  conclude  that 

B  =  V  f  Bg  Bu\  Ur  _  v'  0]  t-7  _  4+ 

b2\  bt, J  l  o  oj 

Note  that  property  (iv)  has  not  been  used  since  it  is  a  consequence  of  (i)-(iii). 

Application  to  Least-Squares  Estimation 

Let  us  now  return  to  the  linear  system  of  Eq.  (27)  and  using  the  notation  defined 
in  Eq.  (4G),  we  are  lead  to  the  system  4c  =  b  where  .4  is  usually  singular.  Hence,  the 
system  may  be  inconsistent,  at  least  numerically,  and  even  if  it  is  consistent,  there  are 
infinitely  many  solutions.  As  suggested  in  the  beginning  of  this  section,  we  will  look  for 
the  (unique)  least-squares  solution  with  minimum  norm:  and  by  this,  we  mean  the  solu¬ 
tion  of  the  extremal  problem  in  Eq.  (48).  The  following  result  gives  the  solution. 

Theorem  3.  Let  A*  be  the  Moore- Penm.-.-  ^.sucJoimvrse  of  A  and  set  c*  =  ATb. 
Then 

||.4c*  -  b||p  —  min  ||.4c  -  b  j  j  2 

C 

and  ||c* ||/>2  =  min { ]•  c [| ^2 :  ||.4c  —  b|| ^2  =  nmi  jj.-lc  —  b j- ^2 } . 


Proof.  Let  .4  =  LEV’7  be  a  singular  value  decompositoin  of  .4  as  described  in  the  sub 
section  entitled  Singular  Value  Decomposition  and  the  Moore- Penrose  Pseudoinverse. 
Then 

||.4c-b|P=||CEV'rc-b|j(, 

=  ||S!'Tc  -  t'7'b||,= 

since  U  is  unitary,  so  that  ()f’a||/-2  =  ([ a [| ^2  for  any  vector  a.  Write 


V'rc  =  b: 


and  UTb  —  {.i] . ju  1  . 


;V(V7C)  -  UTbfn  =  (IT,-.,  -  d,  r  +  •■■  +  (rr,-;4.  -  .h  d 


ir  .4 - +  d-. 

k ■  -*-  1  1  ri 


Hence,  the  minimum  of  ||.4c  —  bjb?  is  attained  at  c  =  <<. 

V1  c  -  . ~  dj  .  'o  4  1 . with  arbitrary 


.  rn  1  where 

.  - 1 .  In  addition,  duce 
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the  solution  c  with  minimum  7 2  norm  is  attained  when  ^*+1  =  •  •  •  —  *  r  =  0.  or  at 
c  =  c*.  where 

VTc*  =  (— di,...,—  &-,0,...,o) 

\<7i  <7*  / 


or  equivalently, 

c*  =  VH+UTb  =  A+b. 

This  completes  the  proof  of  the  theorem. 

Algorithms  for  determining  singular  value  decompositions,  and  hence  Moore-Penrose 
pseudoinverses  A+,  are  available  in  the  literature  [6  -  12],  for  example. 

ESTIMATION  OF  TIME- OF- ARRIVAL 

As  discussed  in  the  subsection  entitled  Spline  Representation  of  Underwater  Acous¬ 
tic  Signals,  the  spline  functions  of  degree  k  defined  on  the  time  interval  [0.  d]  with  knot 
sequence 

t/l-  ^0  ^  <  **•  <C  Uj_  ]  <C  '  t 

where  0  <  to  <  d,  t.n  —  d ,  and  t}  =  tj-\  +  h  for  j  —  1. . . . ,  n  +  k,  with 

h  =  ,  (62) 

n 

will  be  used  to  represent  underwater  acoustic  signals.  As  in  that  subsection,  the  least- 
squares  fit  is  used  to  determine  the  spline  curve.  Since  a  spline  function  is  given  by  the 
spline  series  in  Eq.  (11),  it  vanishes  identically  on  [0,/o]  and  “takes  off’  at  to  (see  Fig. 

9).  It  is  clear  from  this  model  that  the  initial  knot  1 0  represents  the  time-of-arrival. 
Hence,  this  knot  must  also  be  determined.  While  both  7?  and  d  are  fixed,  with  d  de¬ 
noting  the  length  of  the  time  interval  and  n  the  number  of  interior  knots  in  the  time 
interval,  the  relationship  shown  in  Eq.  (62)  implies  that  determining  t„  is  equivalent  to 
determining  h.  Since  the  larger  the  value  n  is  used  the  better  the  estimation  becomes, 
it  is  advisable  to  choose  a  relatively  large  value  of  n,  provided  that  the  computational 
time  is  reasonable.  Hence,  in  the  mathematical  model  in  Eq.  (13).  the  minimization 
must  be  taken  not  only  over  the  spline  coefficients  c0 . r„_i.  but  also  over  the  non¬ 

linear  parameter  h.  Furthermore,  there  are  at  least  two  reasons  that  we  should  choose 
the  minimum  value  of  those  h  that  solve  the  optimization  problem.  Fimv  if  the  "data 
function"  happens  to  be  piecewise  linear  with  equally  spaced  knots,  then  the  minimum 
h,  or  maximum  1 0.  is  the  exact  time-of  arrival,  while  certain  smaller  estimates  of  6)  still 
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reproduce  the  signal  (see  Example  1  in  Appendix  II).  In  addition,  since  h  is  the  distance 
between  two  consecutive  knots,  the  smaller  the  value  of  h.  the  better  the  approximation 
of  Sk  =  Sk,h  [given  by  Eq.  (11)]  to  the  measured  signal  /  is  obtained.  More  precisely,  we 
will  study  the  following  extremal  problem.  Let 

n  -  1  /  n  -  1  \ 

A't(M)  =  ||/  -  .  (03) 

J=»  \j  =  <>  / 

where  c  =  (cq.  . .  .  .  r„-i )  and  A  >  0  is  fixed.  Determine  the  set  H  ~  {h}  of  values  h  such 
that 


7v *(  f>.  c)  =  inf  A\( h.  c)  (G4) 

h  ,c 

for  some  sequences  c  =  (c0, .  .  .  ,  c„_  { ),  where  inf  denotes  the  infimum  ior  "minimum" ) 

and  is  taken  over  all  possible  h  >  0  and  all  c  =  (c„ . ).  If  should  be  emphasized 

that  the  minimization  is  taken  independent  of  the  order  of  h  and  c.Let 

/t*  =  h(k)  =  inf  H 

be  the  greatest  lower  bound  of  the  set  H  =  {h}.  Then  by  E<p  (G2).  the  time- of- arrival  of 
the  acoustic  signal  with  measurement  /  is  given  by 


to  —  d  —  nh * 


(Go) 


Existence,  Uniqueness,  and  Characterization 

As  discussed  in  the  subsection  entitled  Spline  Representation  of  Underwater  Acous¬ 
tic  Signals,  both  the  L 2  and  f2(w)  norms  will  be  used.  Given  any  >  0,  by  the  r:  -frui¬ 
tion  of  infinum,  there  exists  a  pair  (/?o,c0)  where  //()  >  0.  such  that 

A\(//,).c())  <  inf  Kk(h.c)  +  c.  (GG) 

h  ,c 

Since 

I\k{h0.  Co)  >  inf  I\k{h0.c) 

C 

>  inf  {inf  A\.(/t.c)}. 

*  >  0  c 

we  have 

inf  (inf  Kk{h.  c)l  <  inf  A\(  h.  e)  +  £ .  i  G7 ) 

/i>0  1  c  J  h.c 

Since  this  inequality  holds  for  any  f  >  0.  we  may  conclude  that 

inf  ( inf  A'*(  //.  c)l  <  inf  A’ k{h.  r).  (  GS  1 

h  >0  l  c  J  h  ,c 

However,  since  it  is  clear  from  definition  that  the  quantity  on  the  right-hand  -ide  A  no 
greater  than  that  on  the  left-hand  side,  it  follows  that  the  extremal  problem  we  wanted 
to  solve  becomes  an  iterated  extremal  problem,  namely: 
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inf  Kk(h,c)  =  inf  <  inf  A'*(  h.  c)  >  .  (69) 

h,c  h>0  l  c  J 

Now,  for  each  fixed  h  >  0,  since  the  extremal  problem  of  finding  a  c(h)  such  that 

Ek(h)  :=  Kk(h,c(h))  =  inf  Kk(h.c)  (70) 

C 

is  the  (linear)  least-squares  problem  discussed  in  the  subsection  entitled  Spline  Repre¬ 
sentation  of  Underwater  Acoustic  Signals,  which,  as  we  have  seen,  always  has  a  solution 
by  solving  a  system  of  linear  equations,  and  in  fact,  this  solution  is  always  unique  for 
L1 ,  and  also  unique  for  £2(w)  if  the  minimum  norm  (or  Moore-Penrose  pseudoinverse) 
solution  is  used,  we  see  that  the  original  extremal  problem  in  Eq.  (64)  will  also  have  a 
solution  provided  that  there  exists  an  h  >  0  such  that 

Ek(h)  =  inf  Ek(h).  (71) 

hy*  0 

But  the  existence  of  h  is  clear,  since  Ek(h)  is  a  continuous  function  on  0 ,d/n  and 

Ek(  0)  =  !|/||2 

cannot  be  a  minimum.  Let  H  be  the  non-empty  set  of  h  >  0  that  satisfies  Eq.  (71).  The 
greatest  lower  bound  of  H ,  which  is  clearly  a  positive  number  denoted  by  h* ,  is  unique. 
Summarizing  the  above  argument,  we  have  the  following. 

Theorem  4.  There  exist  a  unique  h*  >  0  and  a  sequence  c*  =  (cj, . . . ,  c*  _•, )  such  that 

Kk(h*,  c*)  =  inf  Kk(h,  c). 

h,c 

provided  that  the  minimum-normed  least-squares  solution  is  used  when  the  C2( w)  norm 
is  considered,  where  the  quanity  I\k(h,  c)  is  defined  as  in  Eq.  (63),  and  h*  is  the  greatest 

lower  bound  of  all  h  >  0  that  satisfies  Eq.  (64).  Furthermore,  (h*,c*)  can  be  achieved 
as  follows.  For  each  h  >  0,  let  c(  h )  be  the  unique  (minimum-normed)  least-squares  solu¬ 
tion  of  Eq.  (70).  Find  the  set  H  of  absolute  minima  of  the  continuous  function  Ek(h)  on 
[0,d/n].  (An  absolute  minimum  must  also  be  a  relative  minimum  here.)  Then 

h*  =  inf  H, 

and  c*  =  c*(/i*)  is  the  (minimum-normed)  least-squares  solution  of  Eq.  (70)  for  h  —  h* . 

Recall  that  the  time-of-arrival  is  given  by  to  =  d  —  nh‘  in  Eq.  (65).  The  main  dif¬ 
ficulty  in  *he  procedure  outlined  above  is  solving  the  nonlinear  problem  of  finding  the 
absolute  minima  of  Ek(h ):  namely,  the  extremal  problem  in  Eq.  (71).  For  the  continu¬ 
ous  setting  (or  L2  norm),  this  procedure  will  be  greasy  simplified  by  giving  an  explicit 
formulation  of  the  error  functional  Ek(h). 

Estimation  for  the  Continuous  Setting 

For  any  discrete  da » a  -  1 . A  and  0  =  tq  <  •  •  •  <  rv  =  <1.  let  /(/) 

be  the  piecewise  linear  function  on  [0,d],  linear  on  each  subinterval  [r,.rl  +  |],  such  that 


Chili 


f(r,)  =  /,.  Let  be  the  n-vcctor  as  defined  in  Eq.  (21).  By  using  the  same  change  of 
variables  as  in  the  suhseetion  entitled  Computation  of  the  Coefficient  Matrices,  we  have 


fO 

r/i 


r  f°  i 

Jh,  0 

rO 

.  J  h  —  1  . 


=  h  f  /(  h{  t  -  n  +  / )  -f  d)AT(  t  )<!t. 

Jo 


for  i  —  (), .  .  .  .  a  —  1.  Returning  to  Eqs.  (20),  (30)  and  (33).  let  .4*  /,  =  where 


h 


k 

‘J 


h  /  Nk(t)Xk(t  +  i-j)dt. 
Jo 


(72) 


isj  =(),...  ,n  —  1.  Then  c  =  c (h)  is  the  (unique)  solution  of  the  system  of  linear  equa¬ 
tions: 


To  study  this  system,  set 


(Ak,h  +  A/„  )e  —  f// . 


(73) 


D  = 


ATCiATC  +  >  -  j)<lt 


<)<i.j<n-l 


f(h(t  —  t>  +  /)  +  r/)AT  (t)dt 


V  L-/0  J  0 < i < « —  1 

It  is  clear  that  D  is  a  positive  definite  symmetric  n  x  u  matrix.  Let 


(74) 


A 


'A, 

.0 


/•"N 

o 


A, 


(75) 


he  the  diagonal  matrix  of  eigenvectors  of  D  where  A]  >  .  .  .  >  A„  >  0  and  form  a  unitary 
matrix  P  from  the  corresponding  eigenvectors.  Then  D  —  P  f\  P 1  .  Hence.  Eq.  (73) 
becomes: 


( ftP  A  P1  +  A/„)c  =  //  b 


or 

P(h  A  f  A r„)PTc  =  hb. 

so  that 

c  =  PTP'b 

where 


li\  I  U 


/,  A  f  A 


( 7G ) 


(  7Si 
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Note  that  both  b  =  b (h)  and  c  =  c (h)  are  functions  of  h.  The  function  Ek(  // )  in  Eq. 
(70)  to  be  minimized  now  becomes: 


Ek{  h)  =  hcr Be  4-  Ac1  c  —  2/;cr  b  + 


r  i/( 

Jo 


To  simplify  the  computational  procedure,  it  is  recommended  to  set 

b  =  PT  b 

where  b  is  defined  in  Eq.  (74).  Then  Eq.  (80)  can  be  simplified  to  be 


Ek(h)=  (  \h')\2dt-Y 

,=  1 


h  A,  +  A 


where  b  =  (bo . bn-i)-  It  is  now  clear  that  Ek(h)  is  differentiable  with  continuous 

derivative  E'k(h )  given  by 


i ,  i  \  ST'  h2^,  +  2frA -2 

=  -  E  (hXi  +  A)4- 


V'  h  V 


It  is  also  interesting  to  note  that  E'k(0)  =  0.  Indeed,  h  —  0  gives  a  relative  maxi¬ 
mum  Ek(h).  To  determine  the  objective  function  Ek(h)  in  Eq.  (S2).  we  must  compute 
the  eigenvalues  {A] , . . . ,  A„  }  of  B  =  Bn  and  the  unitary  matrix  P  formed  by  the  cor¬ 
responding  (orthonormal)  eigenvectors  {i/j, . . . ,  u„},  namely:  P  =  [« j  . . .  v„\.  In  the 
following,  we  give  the  example  where  k  =  1. 

Linear  Spline  Estimate 

From  Eq.  (37),  we  see  that  for  k  =  1, 


=  B„  =  f  f 

.  Jo 


l\\(t)Nx(t  +  i  —  j  )dt 


0  <  i  ,j  <  n  —  1 


'4  1 

1  4  1 


1  4  1 

1  2. 

From  Gershgorin's  theorem,  there  is  one  eigenvalue  in  the  interval  [3.5],  u  —  2  eigenval¬ 
ues  in  the  interval  [2,6],  and  one  eigenvalue  in  the  interval  [1.3].  We  can  be  much  more 

precise  by  determining  the  eigenvalue  and  eigenvector  pairs  (A,,  u, )./  =  1 . n.  more 

explicitly.  To  do  so.  let,  us  consider  the  homogeneous  linear  system,  where  we  have  mul¬ 
tiplied  the  matrix  D  by  6: 

■(4-//)  1  I 

1  (4-//)  1  [•"•!  PM 


1  (4-p)  1 

1  ( 2  -  /0 . 


C'lmi 


By  setting  x  =  $  —  2.  Eq.  (S5)  can  he  reformulated  as: 

'  -2  J-iji  +  y2  =  0 

y  i  —  2.11/2  +  i/.i  =  0 


< 


(SC) 


y ii  —  2  2 jyn  —  i  "i"  i/n  0 

.  j/,,-1  -  2(J  +  1  )y„  =  0. 

It  is  well  known  that  the  two  linear  independent  solutions  of  the  second  order  difference 
equations  i)k+>  —  2.V’ik  +  i  +  xJk  =  0  are  the  Chebysliev  polynomials  of  the  first  and  second 
kinds.  Tk(.r)  and  respectively.  where 


T^(x )  =  cos  k  9 

sin(  k  +  \  )9 


I'kir)  = 


sin 


9 


(Si) 


with  9  =  cos-1  x.  0  <  9  <  7T.  Since  C*0(.r)  ~  1  and  L  \(x)  =  2.r  we  have  yx  =  f  r-<(.r). 

A1  =  1 . n.  so  that  all  except  the  last  equation  in  Eq.  (SG)  are  satisfied.  To  satisfy  the 

last  equation  in  Eq.  (S6).  we  must  have: 


or  equivalently, 


sin(n  —  1)0 
sin  6 


(2x  + 


1) 


sin  n  9 
sin  9 


=  0 


(SS) 


sin( n  +  1)0  +  2  sin  u0  =  0  .  0  <  9  <  (S9) 

This  conclusion  implies  that  [,rj  <  1  or  2  <  /x  <  6.  However,  as  pointed  out  above,  it 
follows  from  Gershgorin's  theorem  that  there  is  an  eigenvalue  in  Refs.l  and  3  and  hence, 
it  is  possible  that  1  <  //  <  2.  To  determine  //  with  1  <  //  <  2.  9  may  be  considered  to  be 
complex.  This  is  valid  since  i/*  =  f  T_i(.r)  where  x  =  cos  9.9  complex,  still  satisfies  Eq. 
(SG).  A  careful  investigation  reveals  that  a  complex  0.  which  gives  a  real  x  <  —1.  is  the 
only  root  other  than  the  (n  —  1)  real  roots  of  Eq.  (SO)  in  (0.  r  ). 

Case  1.  Real  Solutions  0  That  Give  the  First  ( n  —  1 )  Positive  Eigenvalues: 


Let 


Then  we  have  : 


f(9)  —  sin( n  +  1  )0  +  2 sin  ii9. 


=  2  sin 


u  +  1 


=  21  -1  )J  +  ]  sin 


I  00) 


so  t hat  for  j  =  1 . ii. 


>  (>  1 1 


01 


2$ 
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/  _  ((n  ~  1)7r  (7l  ~ 

V  n  +  1  n  ) 

do  not  overlap.  Note  also  that  these  intervals  all  lie  in  (0,  ~),  and  in  view  of  Eqs.  (91) 
and  (92),  there  is  one  root  0i  of  Eq.  (89)  in  each  //, /  =  1, - n  —  1. 

Case  2.  Complex  Solution  0  That  Gives  the  nth  Positive  Eigenvalue: 

Let 

g(t )  —  sinh(n  +  1  )t  —  2sinh  nt.  (94) 

Observe  that  since  there  must  be  an  eigenvalue  p  in  (0,2),  or  x  <  —  1.  we  set  6  =  tt  +  jt 
where  j  =  y/—\.  This  gives 

+jt)  =  (-1)"+1  sin^;(n  +  l)fj  -f  2(-l)n  sin (jnt) 


=  (  —  l)"+1j  sinh(n  +  1  )t  —  2sinh7?f 

=  (-i  )n+ljsit)- 

That  is,  we  must  solve  for  the  real  f  in 

()( t )  =  sinh(  ??  +  1  )f  —  2  sinh  nt  =  0, 0  <  t  <  tc  (95) 

and  the  (unique)  root  tn  [or  complex  root  Bn  =  tt  +  jt„  of  Eq.  (S9  7  gives: 

.r „  :=  cos(  ~  +  j6)  =  —  cosh  t„  <  —  1  ( 9G  i 

or 

A„  :=  4  +  2.rn  <2.  (9<  ) 

We  can  now  make  the  following  conclusion.  Let  €  I\.1  =  1.  .  .  .n  —  1,  be  sobuunis 

of  Eq.  (S9)  and  t„  be  the  solution  of  Eq.  (95).  Then  for  I  =  1 . "  —  1.  recalling  the 

factor  1/6  in  Eq.  (84).  the  eigenvalue-eigenvector  pairs  of  the  ma’iix  B  are. 


29 


Chui 


A|  =  (4  +  2  cos  )/G 

1 


U,  = - 

and  for  /  =  n.  the  eigenvalue-eigenvector  pair  is 


i  /_» 


sin  9 1 
sin  2  9i 

L  sin  nBi  j 


'  An  =  (4  -  2  cosh  t„  )/G 

—  sinh  t„ 

1 

sinh  '2t„ 

it  - 

"  t  \1/J 

.i—l  i”  sinh  nt „  . 

I  98 


99 


Hence,  for  k  —  1  (the  linear  spline),  the  following  algorithm  can  he  used  to  estimate  the 
time-of-arrival.  In  all  the  algorithms  presented  in  this  report,  it  should  be  remarked  that 
the  original  extremal  problem  can  lie  written  as  a  nested  extremal  problem,  as  verified 
in  the  subsection  entitled  Existence.  Uniqueness,  and  Characterization.  That  is.  a  linear 
regression  is  first  performed  by  simple  linear  algebra  to  determine  the  spline  coefficients, 
and  then  a  nonlinear  optimization  procedure  follows. 


Algorithm  I  (Linear  Spline  Estimation  of  Time-of-Arrival  Under  Low 
Noise  Condition). 

( 1° )  Compute  A] . A„_| .  A„  and  U| . U.,_  | .  U„  in  Eqs.  (98)  and  (99). 

(2°)  Let  P  =  [Uj  ...  U„]  and  compute  b  =  P1  b  when'  b  is  given  by  Eq.  (74). 

(3°)  For  exact  data,  use  A  =  0.  while  the  noisier  the  data,  the  larger  positive  A  is  re¬ 
quired.  Fix  a  A  >  ()  and  determine  the  set  of  h  such  that 

E\(h )  =  mm/)>o  E,  (h ).  where  E\(h)  is  given  by  Eq.  (82)  with  /i,_j  being  the 
ith  entry  of  b  (which  depends  on  h). 

(4°)  Determine  the  smallest  value  h*  >  0  among  all  values  h. 

The  time-of-arrival  is  given  by  t =  d  —  n  h* .  [See  Eq.  (05)].  The  choice  of  A  as  a 
function  of  the  noise  of  the  signal  will  be  studied  in  the  forthcoming  report  [4], 

In  computing  A[ . A„_|.  the  “bi-section"  method  may  be  used  to  determine  0,  £ 

I,  since  the  values  of  f(d)  have  opposite  signs  at  the  two  end-points  of  I,.  To  compute 
A„.  Newton’s  method  may  be  used  to  search  for  the  unique  root  t„  >  0  of  Eq.  (95). 

In  computing  the  data  vector  b.  and  hence  b(b  =  P1  b).  note  that  it  dejiends  on  h. 
In  fact,  since  k  =  1.  the  {>  +  If  component  of  b  is  given  by 

h,  =  jf  tf(  h(  t  -  7)  +  ?)  -f  d)dt  +  j  (2  -  -  »  +  /)  +  d)dt  (100 

for  7—0 . 77  —  1 .  where 


f  1  fol  7=0 . 77  -  2 

If)  for  7  -  77  -  1 
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The  integrals  in  Eq.  (100)  can  also  be  written  as 

fd—(n  —  i—l)h 


b‘  =  hL.-m 

i  i‘d-(n-i-2)h  / 

+  f>iy  /  (2- 

11  J  d—{n  —  i  —  l)h  \ 


f(t)dt 


i  -  d 


n  +  ^  f(f)dt. 


( 101 ) 


Here,  the  factor  1/h  which  gives  a  factor  of  1  /h2  for  bt  should  not  be  evaluated  since  it 
cancels  with  h 2  in  E^(h).  [See  Eq.  (82)].  Recall  that  f(t)  is  the  piecewise  linear  function 
determined  by  the  interpolation  condition  f{r})  =  f },  where  ( r}.fj )  is  the  data  set. 

In  determining  h  in  (3°),  instead  of  minimizing  the  quantity  E\(h ),  it  is  equivalent 
to  maximizing  the  quantity 


T<'->  =  £r 


i=i 


h\,  +  A6*'1 


(102) 


and  as  pointed  out  above,  the  numerator  h2  is  cancelled  out  with  the  factor  1/h2  from 

/S  /s 

b2_ j.  Recall  that  b ,  depends  on  h.  To  facilitate  the  optimization  process,  the  derivative 
of  T(h)  may  be  used.  Let  c  =  [c0  . . .  en_ j ]  and  c  =  [c0  . . .  c„ _ j )  =  P  c  where 


•  d~  (n  —  i  —  1  )h 


d  —  t 


C;  :  = 


J d  —  ( n  —  i )  h  h 

.d-(n-t -2)h  d_t 


r  (l  — (  n  —  l  —  4 

Si  \ 

J  d  —  (n— i  —  1 ) 


f(t)dt 
f{t)dt 


h 2 


-(n—  i  —  1  )h 

—  ( 1  —  6,  )(n  —  i  —  1 )/( d  —  ( n  —  i  —  1 ) h ). 


(103) 


Then  we  have 


*”<*)  =  £ 


-A, 

(/*A,  +  A)2 


2 

h,\,  T  A 


(104) 


Estimation  with  Splines  of  Arbitrary  Degree 

To  estimate  with  spline  functions  of  higher  degree,  the  method  derived  in  the  sub¬ 
section  entitled  Linear  Spline  Estimate  cannot  be  applied:  and  hence,  we  must  depend 
on  numerical  estimates  of  eigenvalue-eigenvector  pairs  directly.  Let 


D 


k- 

n 


■X  k(t  )A  *(/  +  i 


-  j  )df 

.  IK'  #.;«*'  r,  -  1 


i  i  05  1 
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where  t lie  integrals  are  computed  using  the  procedure  derived  in  the  subsection  entitled 
Computation  of  the  Coefficient  Matrices.  For  instance,  when  cubic  splines  are  used  we 
have 


D 


i 

n 


'C„_, 

D1 


1 10G  j 


where  Cn--\.D.  and  E 3  are  given  in  Eqs.  (41).  (42),  and  (43).  respectively.  Let 


be  the  eigenvalue-eigenvector  pairs  of  77*  where  A,'  >  A*  >  .  .  .  >  A*  >  0  and  each  C/  is 
normalized  to  have  unit  length.  The  algorithm  to  determine  the  time-of-arrival  by  using 
a  k,h  degree  spline  curve  can  be  described  as  follows. 

Algorithm  II  (Estimation  of  Time-of-Arrival  by  Splines  of  Arbitrary 
Degree  k  under  Low  Noise  Condition). 

( 1°)  Choose  k  (depending  on  the  desirable  smoothness)  and  compute  the  Bernstein 
coefficients  of  the  77-spline  Xk  as  in  the  subsection  entitled  Computation  of  77- 
Splines.  (See  Figs.  G  and  7). 

(2°)  By  using  the  formed,*  in  Eq.  (37)  and  following  the  procedure  described  in  the 
subsection  entitle-’  Computation  of  the  Coefficient  Matrices  compute  the  ma¬ 
trix  .4*  [If  „  ibic  spline  curve-  is  used,  skip  these  two  steps  and  use  the  for¬ 
mula  in  F  t.  ,U).j 

(3°)  C'ompub  the  eigenvalue  and  eigenvector  pairs  (Af.U*).t  =  1 . 11,  and  A*  > 

.  .  .  >  a))  >  0.  of  the  matrix  Ak  /, .  (E.g.  the  routines  in  Ref.  13  may  be  used.) 

(4°)  Let 

p=  [uf...u;;] 

and  compute  b  =  Prb  where  b  is  given  by  Eq.  (74). 

(5°)  For  exact  data,  use  A  =  0.  while  a  positive  A  may  be  used  for  noisy  data.  Fix  a 

A  >  0.  and  determine  the  set  of  h  such  that 


Ek[h)  -  min  Ek(h). 

h>  0 

where  Ek(h)  is  given  by  Eq.  (82)  with  being  the  i,h  entry  of  b  (which  de¬ 
pends  on  h ). 

(G°)  Dt-fermine  the  smallest  value  IE  >  0  among  till  values  //. 

Again  a  method  to  determine  the  value  of  A  as  a  function  of  the  noise  will  be  di-- 
nnsed  in  the  forthcoming  report  [4j.  For  the  time  being,  use  A  =  .01  or  even  a  "m;ill'-: 
value  if  the  noise  level  is  very  low. 

Estimation  for  the  Discrete  Setting 

If  the  noise  level  is  fairly  high,  algorithms  I  and  II  are  not  applicable  sine*  the;,  i- 

no  reasonable  criterion  to  determine  the  data  function  /.  In  this  case  the  I2  noim  wit);  a 
relatively  larger  value  of  A  should  be  used.  The  value  of  A  -  1  is  reeommedned.  although 
it  varies  with  the  noise.  See  the  forthcoming  report  [4]  for  a  better  choice  of  V 
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Let  ( T(,fe).C  =  1, . . . ,  N,  be  a  set  of  data  information.  Set 


Akh=  b 


where  6*  is  given  by  Eq.  (31)  with 


=  Nk  f  n  -  i  -  -(d  -  t)  1  . 


Define 


fh  :  = 


and  let 


£”=  1  f(Bk,th.0(Tf)U'( 


.£"=1  ftBk%th.n-\{Tt)wt_ 


A  —  Akh  +  A/n 


(109) 


(110) 


as  in  Eq.  (46).  Then  we  have  the  following  algorithm  to  determine  the  time-of-arrival  by 
using  a  kth  degree  spline  curve  and  the  72(w)  norm. 

Algorithm  III  (Estimation  of  Time-of-Arrival  with  the  £2  Norm  for  Noisy 
Data). 

(1°)  Choose  k  (depending  on  the  desirable  smoothness)  and  compute  the  Bernstein 
coefficients  of  the  R-spline  Nk  as  in  the  subsection  entitled  Computation  of  B- 
Splines.  (See  Figs.  6  and  7). 

(2°)  Compute  the  polynomial  pieces  of  Nk  by  using  the  formula  in  Eq.  (5)  with  a*) 
being  the  Bernstein  coefficients  and  (pkj(u,v)  defined  by  Eq.  (4),  where  u,  v  are 

given  in  Eq.  (2)  with  [a,  6]  being  the  corresponding  interval  of  the  polynomial 
piece.  Set 


=  Nk  (n-i-  d  -  f)^ 


(111) 


(3°)  Compute  bk}  in  Eq.  (31)  and  /,  by  using  uy  =  1.7  =  1, . . . ,  jV,  and  0  <  i,j  < 
n  —  1. 

(4°)  Fix  a  positive  value  A,  say  A  =  1  in  Eq.  (110).  (The  noisier  the  data,  the  larger 
the  value  of  A  is  recommended.)  Determine  an  SA  D  (singular  value  decomposi¬ 
tion)  A  —  UT,VT  of  .4,  where  U  and  V'  are  unitary  and 


where  rr.  >  .  .  .  >  rr„,  >  0. 
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(5°) 


Let  u,  ho  the  ith  column  of  U  and  v,  the  i"‘  column  of  U. 
[c(*  .  .  .  r*_j]  by  using  the  formula 


Compute  c*  = 


c 


* 


(113) 


where  f/,  = 
(6°)  Compute 


/o  •  ••  fn-l 


has  been  computed  in  (3°). 


i«h)  =  Y. 

i=i 


(114) 


Here,  we  have  used  ic,  =  1.  The  value  of  A  must  be  the  same  as  the  A  in  (4°). 
(A  smaller  or  larger  value  of  A  is  used  depending  on  the  noise  level  of  the  data.) 

(7°)  Determine  the  set  of  h  such  that 


K  (7i)  =  min  K(h ). 
V  /  /i>o 


(115) 


(8°)  Determine  the  smallest  value  h*  >  0  among  all  values  h. 

(9°)  Compute  t0  =  d  —  nh* . 

Then  to  is  the  time-of-arrival.  Note  that  algorithms  from  Refs.  C  to  9  can  be  used 
in  step  (4° ). 
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